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Abstract 

We report here on an extension of a previous study Kirsh et al. (2009) of planetesimal- 
driven migration using our iV-body code SyMBA (Duncan et al., 1998). The previ- 
ous work focused on the case of a single planet of mass M em , immersed in a plan- 
etesimal disk with a power-law surface density distribution and Rayleigh distributed 
eccentricities and inclinations. Typically 10 4 to 10 5 equal-mass planetesimals were 
used, where the gravitational force (and the back-reaction) on each planetesimal 
by the Sun and planet were included, while planetesimal-planetesimal interactions 
were neglected. The runs reported on here incorporate the dynamical effects of a 
gas disk, where the Adachi et al. (1976) prescription of aerodynamic gas drag is 
implemented for all bodies. In some cases the Papaloizou and Larwood (2000) pre- 
scription of Type - I migration for the planet are implemented, as well as a mass 
distribution. 

In the gas-free cases, rapid planet migration was observed - at a rate independent 
of the planet's mass - provided the planet's mass was not large compared to the 
mass in planetesimals capable of entering its Hill sphere. In such cases, both inward 
and outward migrations can be self-sustaining, but there is a strong propensity 
for inward migration. When a gas disk is present, aerodynamic drag can substan- 
tially modify the dynamics of scattered planetesimals. For sufficiently large or small 
mono-dispersed planetesimals, the planet typically migrates inward. However, for 
a range of plausible planetesimal sizes {i.e. 0.5 - 5.0 km at 5.0 AU in a minimum 
mass Hayashi disk) outward migration is usually triggered, often accompanied by 
substantial planetary mass accretion. The origins of this behaviour are explained 
in terms of a toy model. The effects of including a size distribution and torques 
associated with Type - I migration are also discussed. 
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Planetary Dynamics 



1 Introduction 

Currently, over 460 extrasolar planets are knowrPI. along with over 40 sys- 
tems containing multiple planets. Most of the extrasolar planets detected to 
date have masses co mparable to that of Neptune, or larger. Furthermore in a 



recent summary by lUdry et al.l (120071 ). at least ~ 6% of stars surveyed have 



giant planets interior to ~ 5.0 AU, so giant planets appea r fairly common in 



stellar systems. Moreover, results from the HARPS survey ISousa et al.l (120081 ) 
shows that Neptune-mass extrasolar planets are found in ~ 40% of the stars 
surveyed. Most likely, these giant planets formed via a similar process that 
formed the four giant planets in the Solar System. 

However, the masses and orbits of these extrasolar planets display a wide vari- 
ety of configurations: e.g. Neptune and Jupiter-mass planets with short orbital 
periods, isolated planets with large orbital eccentricities, multiple planet sys- 
tems in resonance, and planets orbiting components of stellar binaries. Several 
analytical models have been proposed to explain the various aspects of planet 
formation, but most of these have not been tested numerically. Until recently, 
very littl e had been done on gia n t planet core format ion using iV-body sim- 
ulations ( Thommes et al. . 2003 ). Levison et al. ( 2010f ) (hereafter referred to 



as LTD10) recently completed a comprehensive set of computer simulations 
which included a number of physical processes that might enhance accretion 
onto planetary embryos. As discussed in Section 2, the most successful models 
were those in which one or more embryos spontaneously underwent a burst of 
outward migration induced by planetesimal scattering. 

In an attempt to further our understanding of some of the results in LTD10, we 
are undertaking a detailed investigation of the combined effects of planetesimal 
scattering and aerodynamic drag on the growth and evolution of giant planet 
cores. Our goal in this paper is to understand the case of the dynamics of a 
single core interacting with a disk of planetesimals and gas. In what follows, 
we provide some background in §2, then briefly discuss our implementation of 
the relevant forces in §3. In §4 we discuss the results of simulations including 
aerodynamic gas drag, for a disk of mono-dispersed planetesimals. A toy model 
which explains the results is presented in §5. The effects of a planetesimal size 
distribution is presented in §6, and Type - I migration in §7. A summary and 
conclusion is presented in §8. 



See http://exoplanet.eu 
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2 Giant Planet Formation 



The formation of giant planets in the widely adopted core accretion model 
can typically be described in four stages. The first stage involves the forma- 
tion of planetesimals, which we do not model in this study. The next stage 
involves the runaway accretion of planetesimals by a small fraction of those 
planetesimals which happen to grow a bit larger, and then grow much faster 
than all the others (jWetherill and Stewartl . Il989[ ). When these large bodies 
become sufficiently massive and well-spaced such that each dominates the vis- 
cous stirring in its feeding zone, the runaway growth gives way to the oligarchic 
growth stage. An embryo's feeding zone is the annulus about its orbit where 
small bodies can suffer strong gravitational impulses. Typically this feeding 
zone extends from 1.0 - 3.5 Hill radii on either side of the embryo's orbit, 
and is the source of most of the material which the embryo accretes. During 
the oligarchic stage, the large embry os grow in lockstep, maintaining similar 



mass es and uniformly spaced orbits ( IKokubo and Ida . Il998l ; iThommes et al 
20031 ). The final stage in the outer solar system is characterized by the rapid 



accumulation of a gaseous envelope by the embryos; in the inner region it is 
characterized by the giant impact phase of terrestrial planet formation. 



However, the core accretion model has its weaknesses. In particular, the accre- 
tion of a massive atmosphe re requires a solid core of mass ~ 1 M ffl to trigger a 



rapid gas accretion phase (iMizunol . Il980l ; iPollack et all Il996l ; iHubickyj et al. 



20051 ). The difficulties of reaching this threshold are threefold: 



(1) Accretion has to be sufficiently efficient to concentrate enough mass into 
at least one body, and potentially multipl e bodies. 

(2) Accretion has to occur within ~ lOMyr (IHaisch et all 1200 ll ). such that 
there is ~ 10 2 M e left in the nearby disk to furnish an envelope. 

(3) Migration due to embryo-disk tidal interactions (cf. §3.2.2), threatens to 
deposit core-sized bodies into the central star faster t han t hey can accrete 
flWardl . Il986l : iKorvcanskv and Pollack! . Il993l : IWardl . Il997h . 



Several analytical models have been proposed to mitigate these problems, and 
some of these have been tested numerically by LTD10. In particular, LTD10 
numerically integrated the orbits of a number of planetary embryos embedded 
in a swarm of planetesimals. Their simulations included simplified models of 
various combinations of the following effects: (1) aerodynamic drag on small 
bodies, (2) collisional dam ping, (3) extended atmospheres around the embryos 
( Inaba and Ikomal 12003 ). (4) embryo eccentricity damping due to gravita- 
tional interaction with the gas disk, (5) fragmentation o f the planetesimals 
and (6) evaporation and re-condensation at the snow line (ICuzzi and Zahnld . 
20041 ). They found that the gravitational interaction between the embryos and 
the planetesimals generally led to regions near the embryos being cleared of 
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planetesimals before much accretion onto the embryos could occur. However, 
the most successful phases of embryo growth occurred when the gravitational 
scattering of the planetesimals, together with the effects of aerodynamic gas 
drag led to the rapid outward migration of one or more embryos. We show in 
this paper that many of the main features of the embryo-planetesimal inter- 
actions that lead to rapid outward migration and planet growth are demon- 
strated by the single embryo case which we discuss next. 



3 Physical Processes in Circumstellar Disks 

There are several physical processes that can occur in circumstellar disks; 
some of these are only relevant to planetesimals, while others only to larger 
embryo-sized bodies. Specifically, the dynamics of planetesimals and embryos 
will be affected by the gravitational perturbations from other massive bodies, 
as well as gas effects. Radiative forces are not very important for 0.01 - 100 
km size bodies over the timescale under consideration (i.e. ~ lOMyr), so we 
neglect such forces in the subsequent discussion. 

3. 1 Gravitational Effects 

The dominant gravitational influence in the circumstellar environment, for 
planetesimals and embryos, is the central star. However, in the vicinity of 
other massive bodies (e.g. embryos), the gravitational tidal influence of those 
massive bodies will dominate. The transition is characterized by a length scale 
called the Hill radius, which defines a sphere about each body where its grav- 
itational tide dominates the gravitational influence from the central star: 



where M and a are the mass and the semi- major axis of an orbiting body, 
while M* is the mass of the central star. At 1.0 AU an Earth-mass object 
would have Rh c± 0.01 AU, while at 10.0 AU a Jupiter-mass object would have 



In the event of a close encounter, the embryo will tend to scatter the planetes- 
imal to a smaller or larger orbit, exchanging energy and angular momentum. 
Consequently, the embryo will respond by moving in the opposite direction 
of the planetesimal, albeit by a much smaller amount. Since an embryo is 
surrounded by a swarm of planetesimals, it will scatter numerous planetesi- 
mals as it moves along its orbit. Furthermore, if the probability of scattering 




(1) 



R h ~ 0.7 AU. 
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a planetesimal inwards were the same as scattering outwards, there will be no 
net change of the embryo's orbit. However, since the timescale for a scattering 
encounter is slightly shorter inside the planet's orbit, it will preferentially scat- 
ter planetesimals from inside its orbit to outside its orbit. Consequently, the 
embryo will experience a net inward drift, and this inward migration will con- 



tinue so long as t here is sufficient material for it to scatter ( [Fernandez and Ip 



1984 ; iMalhotral. 11993c iGomes et all 20041). This migration is studied in de- 



tail by iKirshl ( 120071 ) and iKirsh et al.l ( 120091 ) in gas-free disks, and we briefly 
summarize their work here. 



In their study IKirsh et al.l ( 120091 ) noted that if a swarm of planetesimals were 
scattered by a much more massive embryo, it could lead to a net exchange 
of angular momentum that would induce the embryo to migrate. The rate an 
embryo's orbital d istance will drift due to planetesimal scattering is given by 



(iKirsh et all 120091 ): 



a 
a 

where P orD is the embryo's orbital period at a em , while M disk = S S oiid( a em) 7r Og m 
is the local mass of the disk, where £soiid(a em ) is the local surface density of 
the solid material in the disk and M & is the solar mass. This rate will be 
independent of M em provided M em <C M enc where M enc ~ 5^Mdi S k is the 
mass in the embryo's encounter region, and £/, = Rh/a cm = (M em /3M^) 1//3 = 
10 _2 (M em /M ffi ) 1 / 3 is the Hill factor. This rate is valid provided that the plan- 
etesimal eccentricities e < 3^, which is always the case considered in this 
paper. 



M, 



disk 
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orb 



Mp 



1 
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3.2 Gas Effects 



3.2.1 Aerodynamic Gas Drag 

While the gaseous component of the disk is still present, solid particles will 
be affected by an aerodynamic gas drag. The strength of this gas drag accel- 
eration is inversely proportional to a particle's radius, so larger particles will 
be less effected. Since embryos can have radii greater than a few hundred km 
(depending on their density), gas drag is not an important effect. 

The dynamics induced by aerodynamic gas drag on a particle is described by: 



_ V ~ Vgas ( v 

^aero — \y) 
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where v and v gas are the velocity vectors of the plan etesimal and the gas , and 
T aero is the aerodynamic gas drag timescale given by lAdachi et al.l (119761 ): 



where p p i sm i and -Rdrag are the mass density and radius of the planetesimal, 
while p gas is the volume gas density and v^ ev is the Keplerian orbital velocity. 
The drag coefficient Cp is a quantity often of order unity, and in general is a 
nonlinear function of the particle's radius and its relative velocity to the gas. 
We discuss the behaviour of Cd in more detail below. 

For the analytical discussion in the subsequent sections, it is useful to define 
damping timescales for planetesimals with small eccentricity e and inclination 
/. Damping rates for e and I of the plane tesimals due to aerodynamic gas 
drag were calculated by lAdachi et al.l (119761 ) : 




-e 2 + -P 
2 



(5a) 
(5b) 



where rj 



( ^gaa V 
V V kep / 



defines the deviation of the gas velocity from the Ke- 



plerian velocity at a given location in the disk. Assuming a power-law volume 



density profile [i.e. p g 



(a) oc a 



and a power-law gas temperature profile 



[i.e. T gas (a) oc a 



l ' 2 } in the disk gives: 77(a) = 6.0 x 1(T 4 (a + ^(a/AU) 1 / 2 . 
The corresponding damping rate for a planetesimal's orbit is given by: 



lrj* + ^ + \P (v + e 2 + ll 2 



(6) 



We now digress temporarily to discuss two quantities that appear in the gas 
drag formulae which are often overlooked. The first is the quantity rj 2 , which as 
it appears in Eq. 5 and Eq. 6, is often significantly smaller compared to values 
of e 2 and J 2 . As a consequence, many studies that include these aerodynamic 
gas drag formulae often employ versions of these equations in the limit when 
i] 2 <C e 2 ,I 2 . However as the aerodynamic drag damps the random motion of 
the planetesimals particularly for small embryo masses and at larger distances 
in the disk, the value of rj 2 will become comp arable to e 2 and J 2 . It seems only 
prudent to implement the full equations from lAdachi et al.l (119761 ). rather than 
their typical approximations. 
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Next is the quantity Co, which is often ass umed to be o f order unit y, and is 
often fixed throughout the disk. However, lAdachi et al.l (119761 ) and iRafikov 
( 120041 ) clearly showed that the value of Cp is a nonlinear function of the 
planetesimal radius and its velocity relative to the gas. Furthermore, the value 
of Cd not only varies over a couple orders of magnitude across several orders 
of magnitude of planetesimal radius, but it also varies with distance in the 
disk. This is illustrated in Fig. 1 which plots the value of Co as a function of 
-Rdrag at 5.0 AU and 10.0 AU, and assumes that a gas volume density varies 
radially as a power-law: p gas = 1.4 x 10 _9 g cm -3 (q/A U)~ n / 4 as pres cribed 
by the Minimum Mass Solar Nebula (MMSN) model (iHarashl [l98lh The 
nonlinear trea t ment of Cd is incorporated in the work of iRafikovl (120041 ) and 



Brasser et aJj (120071). and for t he purposes of our study we implement the 



formulae in iBrasser et al.l (120071 ). 




0.1 



0.01 




'drag 



(km) 



Fig. 1. Behaviour of the drag coefficient Co as a function of planetesimal radius 
-Rdragi assuming a gas volume density of 1.4 x 1CP 9 g cm~ 3 at 1.0 AU which varies 
radially as a power-law: p gas oc a -11 / 4 . It is also assumed that the planetesimal 
density /5 p ismi = 0.5 g cm -3 , which is the same in all the figures. The two different 
curves correspond to two locations in the disk: 5.0 AU (solid line) and 10.0 AU 
(dashed line). 



3.2.2 Type - / Migration 

For embryos, a more important dynamical effect than aerodynamic gas drag 
is the gravitational interaction with the gas disk itself. Embryos embedded 
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in a gas disk will launch density waves at their inner and outer Lindblad 
resonances (which are analogous to mean motion resonances in disks), and as 
a result experience a positive and negative torque from the inner and oute r 
portions of the density wave, respectively ( Goldreich and Tremaine . Il980l ). 
Since the outer torque dominates in locally isothermal disks, this cau ses the 
embr yo's orbit to decay, which is usually termed Type - I migration ( Ward! 
19971 ). The migration rate increases with embryo mass until the embryo is 
large enough to clear a gap in the gas disk, at which point it becomes locked 
into the slower viscous evolution of the disk (i.e. Type - II migration, cf. [Ward 



( 119971 )). However, since embryos do not form a gap until they reach a mass of 
10 - 100 M®, then Type - I migration will dominate embryo evolution until 
the embryo becomes a gas giant core. 



Recent work has revealed that if the gas is not isothermal, then the drag from 
the gas in the horseshoe region can lead to outward driven Type - I migration 



( IPaardekooper and Papaloizoul . 120081 ; iPaardekooper et all I201CH ) . However for 



the radial density (i.e. p gas (cO oc a a ) and temperature (i. e. T ga , R (a) oc a 1//2 ) 



profil es explored in this study, the simple torque formula (cf . IPaardekooper et al 



2 1 Ol . Eq. 47) indicates t hat the embryo's preferred migration direction would 
still be inward. Work by lLyra et al.l ( 120101 ) extends the analysis by including a 
viscously and radiatively evolving disk, which is an important improvement as 
most disk models are assumed static or isothermal. While their simulations do 
indeed show regimes of outward migration, these regimes are typically circum- 
scribed between equilibrium radii (e.g. locations in the disk where the torque 
is zero) which drift inwards as the disk evolves. They find that smaller mass 
embryos are able to decouple from the evolution of these equilibrium radii, but 
since they do not incorporate mass accretion for the embryo this may have 
limited application to our study. 



For Type - I migration, iPapaloizou and Larwoodl (120001 ) give the acceleration 
on an embryo due to the tidal effects 



v v • r v ■ k 

Vtidal = 2- 2 k (7) 

''"a^ypel ' ^e,typel ''"ijtypel 

where r and v are the position and velocity vectors of the embryo, and k 
is the unit vector in the vertical direction, while r atype i, T e ,ty P ei and T ItypcI 
are respectively the decay timescales for the orbital distance, eccentricity and 
inclination: 
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1 + 



ea V 



(9) 



2vrc e \ aj \ M Q J \ M Q 

where S gas = y/7iz s p gas is the local gas surface density, and z s is the scale height 
of the disk. The coefficients c a and c e control the strength of the radial drift and 
circularization of an embryo's orbit, respectively. If the inclination damping 
timescale is not significantly sho rter than the eccentricity damping timescale as 
Papaloizou and Larwoodl (120001 ) argue, then it will not contribute significantly 
to reaching the equilibrium state; thus for simplicity we assume T IjtypeI = 
T e ,typei- Most agree that the choice c e = 1.0 is reasonable; however there is less 
consensus about the exact value for the coefficient c a though there is some 
agreement that c a < 1.0 for the case of an isot hermal equation of state for the 
gas. However it has recently been shown by (IPaardekooper and Papaloizou 



20091 ) that the gas horseshoe torque can substantially modify the timescales 
in Eq. 8, and even reverse the direction of migration, when the gas cooling 
timescale is greater than the horseshoe libration timescale. However we shall 
limit our study to the "classical" Type - I formulae, and leave the inclusion 
of the horseshoe torque for a future study. 



While we can include both aerodynamic and tidal gas interaction prescriptions 
in our iV-body code, we first discuss our simulations where we only consider 
aerodynamic gas drag in the next section. We discuss the inclusion of Type - I 
migration in §7, which will be important for larger embryo masses (e.g. M em > 
1.0M ffi ). 



4 Simulations 



We p erform our numerical integrations using the SyMBA integrator (IDuncan et al. 



19981 ). part of the SWIFT suite of packages. Sy MBA is a mixed- variab l e sym - 



plectic integrator based on the iV-body map of IWisdom and Holmanl (119911 ). 
which has been improved to accurately handle close encounters between par- 
ticles. 



In our simulations we consider two types of particles: massive bodies (e.g. em- 
bryos) which mutually gravitate and can accrete other bodies, and less massive 
bodies (e.g. planetesimals) which do not mutually gravitate or accrete. Using 
a particle-mesh gravity solver, we are able to include the self-gravity of the 
planetesimals. However for the disk masses considered in our simulations, the 
contribution from the planetesimal self-gravity is negligible at best, so we omit 
these calculations in all our simulations. We have implemented the gas drag 
prescriptions as discussed in §3.2, where we assume aerodynamic gas drag and 
Type - I migration applies to planetesimals and embryos, respectively. Since 
we will be concentrating on embryos initially smaller than 1.0 M e , we have 
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omitted the effects of Type - I migration on the embryo for now. We investi- 
gate the effects of Type - I migration in §7, particularly in conjunction with 
the planetesimal - driven migration. 

As a fiducial model, we assumed that the circumstellar disk corresponds to a 
power-law, where the volume density distribution of gaseous material: 



/ a 

VAuy 



V 



In a MMSN, p gas , 



Pgas(O-) 2?) fgPgas,0 

1.4 x 10" 9 g cm -3 at 1.0 AU and a 



(10) 

11/4 flHavashil . 



19811 ). while the disk scale height z s (a) is given by a power-law: 



z,{a) = z sfl i^—j (11) 

where z s ^ = 0.047 AU and f g is a scaling factor of gaseous material in the 
disk, with f g = 1.0 in a MMSN. 

The surface density distribution of solid material in a MMSN is also given by 
a power-law: 



Ssolid(fl) = /s,/snowSsolid,0 ( ^JJ J V"^) 

where £ S oiid,o = 7.0 g cm -2 at 1.0 AU and (3 = a — 5/4 assuming a gas to solid 
material ratio that does not vary radially, while 



I 1.0 ifa<a snow 

J snow — \ .. [i-O) 

(^4.2 if a > a snow 

accounts for the enhancement of material beyond the snow-line, where the 
temperature of the circumstellar disk is cool enough for molecules (e.g. H 2 0, 
CO?, CH4, etc.) to condense into solids. Typically a snow = 2.7 AU (L+fL,^ 1 / 2 







( IChiang and Goldreichl . 119971 ). assuming a condensation temperature of 170 
K. The coefficient f s is a scaling factor of solid material in the disk, with 
f s = 1.0 in a MMSN. 

Each of our simulations involve a single embryo on a circular orbit in the 
mid-plane of the disk, embedded in the centre of a 14Rh wid e annulus of plan- 



etesim als. The planetesimals are given Rayleigh distributed (llda and Makino 



19921 ) eccentricities and inclinations with dispersions a e = 2<jj = 0.004, and 
we assume the density of each planetesimal is p p i sm i = 0.5 g cm -3 . All these 
simulations use 3.2 x 10 4 equal-mass planetesimals, where the planetesimal 
mass is dependent on the disk model and the location of the embryo. 
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In most of our simulations we assume all of the planetesimals have the same 
radius -Rdrag; which is an independent parameter and does not depend on 
the mass or assumed density of the planetesimal. Our planetesimals are ac- 
tually tracer particles, each representing a large number of planetesimals of 
size -Rdrag with similar orbits. However we realize that using mono-dispersed 
planetesimals is not realistic, since we expect mutual collisions among the 
planetesimals to produce a size distribution. It is not feasible to implement a 
self-consistent time dependent planetesimal size distribution at present, but 
we can implement a static planetesimal size distribution: dN/dR^ng oc R^ g 
with 10 _1 km < -Rdrag ^ 10 2 km and in §6 we investigate the behaviour for a 
range of values: 2.0 < 7 < 5.0. 

4-1 Effects of Aerodynamic Drag 

We perform simulations for several combinations of parameters, namely: the 
mass of the embryo M em , the orbital distance of the embryo a em , the radius 
of the planetesimals -Rdrag, the gas volume density f g and the solid surface 
density f s relative to the fiducial MMSN model. In particular we explored 
the following parameter values: M cm = {0.25, 1.0} M e ; a cm = {5.0, 10.0} AU; 
f g = {0.5, 1.0,2.0}; f s = {0.5, 1.0,2.0}; and log {R dTag /km) from -2.0 to 3.0 in 
increments of 0.5. For the two embryo locations at 5.0 AU and 10.0 AU, we 
integrate our simulations for 2 x 10 4 year and 4 x 10 4 year, using time steps of 
0.5 year and 1.25 year, respectively. The results of these simulations are shown 
in Fig. 3 and Fig. 4, which shows the rate of change of the embryo's orbital 
distance as a function of planetesimal radius. The behaviour of the migration 
rate versus planetesimal radius is complicated, though we note four distinct 
migration regimes which are described in Fig. 2 for the case of a M cm = 1.0 M e 
embryo at a cm = 10.0 AU. 

4-1-1 i?drag < 0.01 km 

In this regime, the aerodynamic gas drag is so strong that it causes the semi- 
major axes and eccentricities of the planetesimals to decay very rapidly. In 
fact, almost all the planetesimals tend to stream past the embryo before it 
has an opportunity to interact with them gravitationally. The net result is the 
embryo does not migrate significantly at all. 

4.1.2 .Rdrag ~ 0.1 km 

In this regime, the aerodynamic gas drag is still strong, but the embryo is now 
able to gravitationally scatter some of the planetesimals. Most of the plan- 
etesimals tend to become trapped in an exterior resonance with the embryo, 



11 




Fig. 2. Characteristic orbital evolution of a 1.0 embryo started at 10.0 AU in 
a fiducial MMSN disk (i.e. a = 11/4 and f g = f s = 1-0) for different planetesimal 
radii -Rdrag- Two different simulations are plotted for each -Rdrag> where the only 
difference is that the initial distribution of planetesimals were randomized. 

but the aerodynamic gas drag is constantly damping the random motion of 
the planetesimals and causes their inward migration. Since sufficient mass in 
planetesimals becomes trapped they are collectively able to push back against 
the embryo, and the net result is the embryo migrates inwards along with 
these trapped planetesimals. Such behaviour was seen in the simulations of 
LTD10. 



We can estimate the planetesimal radius that will be trapped in an exterior 
re sonance, an d the r ate at which the embryo will migrate, based on the work 
of IfCary et al.l ( 119931 ) . They quantify the trapping condition through the drag 
parameter, K: 



K = 3PsaS p D (14) 

where they denned the critical drag parameter K tr&v to be where 50% of the 
planetesimals become trapped in an exterior resonance, and ~ 50% stream 
past the embryo. They deter mined the value o f Ktrap numerically for several 
different planetary mass es fcf.lKarv et al.l. ll993l . Table-Ill). If we fit the values 
found in this table from lKary et al.l ~f !993 ). we obtain: 
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Fig. 3. Migration rate as a function of the planetesimal radius -Rdrag f° r an 0.25 
embryo at 10.0 AU in a fiducial MMSN gas disk (i.e. f g = 1.0). The squares, 
circles and triangles respond to simulations to disks with f s = {0.5, 1.0,2.0}, while 
the dashed, dot-dashed and dotted lines refer to the fiducial migration rate due 
to planetesimal scattering in Eq. 2. The migration rate, and its uncertainty, are 
measured from two realizations for each -Rdrag- I n these simulations, the value of Cd 
is computed based on the dynamical and physical properties of each planetesimal. 

(M \°' 73 

WM em ) * 1.9 x 10- 2 AU" 1 ( (15) 
Using this relation we can write -Rdra g ,tra P for a MMSN disk model: 



-0.73 , _ n/4 



We test the predicted values in Eq. 16 of i?drag,tra P by running simulations 
for the case of a 0.25 M e embryo at four locations in a MMSN disk, using 
two choices of Co- fixed to 0.5 and computed based on the physical and dy- 
namical properties of a typical planetesimal. In the latter case, the value of 
-Rdrag.trap in Eq. 16 is solved iteratively since the value of can depend on 
-Rdrag in a nonlinear fashion. To facilitate the determination of the fraction 
of trapped pianetesimals, we assume that the planetesimals are massless in 
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Fig. 4. Migration rate as a function of the planetesimal radius i?drag f° r a 1.0 
embryo at 10.0 AU in a fiducial MMSN gas disk (i.e. f g = 1.0). The squares, 
circles and triangles respond to simulations to disks with f s = {0.5, 1.0,2.0}, while 
the dashed, dot-dashed and dotted lines refer to the fiducial migration rate due 
to planetesimal scattering in Eq. 2. The migration rate, and its uncertainty, are 
measured from two realizations for each -Rdrag- I n these simulations, the value of Cd 
is computed based on the dynamical and physical properties of each planetesimal. 



these simulations. Illustrated in Fig. 5 are the predicted values of -Rdra g ,tra P as 
a function of distance in the disk, for both choices of Cd- Also plotted are 
results of the numerical experiments to verify -Rdrag.trap, and the points agree 
with the predicted values within a factor of three. This discrepancy is more 
pronounced at larger distances, which may be a cons equence of e x trapo lating 
-Rdrag,tra P beyond 5.0 AU where all the simulations of IfCary et al.l (119931 ) were 
performed. 



To estimate the rate at which the trapped planetesimals will push the embryo 
inwards, we compute the value of d aero from Eq. 6 for mono-dispersed plan- 
etesimals of size -Rdrag,tra P , assuming some RMS eccentricity e rms for the plan- 
etesimals. However, the rate at which the ensemble of trapped planetesimals 
and the embryo will migrate inward will be reduced by a factor 1 + M em /M trap , 
where M trap is the mass of the planetesimals that remain trapped in an exterior 
resonance with the embryo. 
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Fig. 5. The critical radius above which a planetesimal will become trapped in an 
exterior resonance with an embryo, as a function of the embryo's semi-major axis. 
The transition radius is compute d from the critical drag parameter K determined 
empirically by Karv et al. (1993), assuming either Cd = 0.5 (dashed line) or Cd 
computed based on the dynamical and physical properties of a typical planetesimal 
(solid line). The points are results of simulations at different locations in the disk, 
squares for Cd = 0.5 and circles for computed Cd- In all the simulations the value 
-K'trap = 6.9 x 10~ 3 AU _1 for a 0.25 Af^ embryo is used, and the planetesimals are 
assumed to be massless. 



^aero(^5 ^rmsi -^drag,trapy 

^ , ro ' trap ~ 1 + M em /M trap 



(17) 



For comparison with Eq. 17 we select the migration rates in Fig. 3, at or 
near -Rdra g ,tra P for two locations in a MMSN disk. From these selected simu- 
lations, we estimate the typical mass of trapped planetesimals M trap for use 
in Eq. 17. Inspection of these simulations reveals that ~ 1-2 M em of plan- 
etesimals remain trapped in the exterior resonance. Plotted in Fig. 6 is the 
measured migration rate of a 0.25 M e embryo at two locations in a MMSN 
disk, for planetesimal populations at or near -Rdrag.trap- Also plotted are the 
value of a ae ro,trap from Eq. 17, where the lines correspond to the assumed RMS 
eccentricity e rms for the planetesimals. We note the good agreement between 
the measured and predicted migration rates for the range of e rms assumed. 
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Fig. 6. The predicted migration rate d acro ,trap of a 0.25 M e embryo as a function 
of the embryo's semi-major axis, that is being pushed by a population of trapped 
planetesimals of size i?drag,trap- This calculation assumes that the planetesimal pop- 
ulation trapped in an exterior resonance contains ~ 0.25 M®, which is consistent 
with the simulations. Plotted are curves of d aer o,trap for different assumed RMS ec- 
centricity e rms of the planetesimal population, in units of the Hill eccentricity e^. 
The circles and squares correspond to the measured migration rate of the embryo 
from our simulations at two different locations, for planetesimal populations at or 

near -Rdrag,trap- 



4-1.3 i? drag ~ 1.0 km 

In this regime, the aerodynamic gas drag damps the random motions of plan- 
etesimals more slowly than that in §4.1.2. While the gas drag is diminished, 
it is still able to shift the relative populations of the scattering material on 
either side of the embryo. This shift in the scattering populations is able to 
tip the imbalance in scattering events from outward to inward, and the net 
effect is the embryo tends to migrate outward in response. This will be dis- 
cussed in greater detail in §5. The speed of the outward migration of the 
embryo is also comparable to the inward migration in the gas-free case, which 
seems to indicate that planetesimal - driven migration is the driving mech- 
anism in the outward migration case. The maximum outward migration rate 
for i?drag ~ 1.0 km appears to be reproduced in most of the combinations of 
the parameters that we examined, with minor variations. 
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4-1-4 #drag > 10.0 km 



In this regime, the aerodynamic gas drag becomes increasingly irrelevant. This 
is evident as the migration rate asymptotically approaches the gas-free migra- 
tion rate, which i n mos t cases agrees favourably with the estimate in Eq. 2 
from iKirsh et al.l ( 120091 ). These are plotted as dashed, dot-dashed and dotted 
horizontal lines in Fig. 4 and Fig. 3. However, the results for the 0.25 M e sim- 
ulations for f a = 2.0 do not agree as favourably with Eq. 2. Since the other 
sets of simulations do ag ree, it is worth noting that the migration rate in Eq. 2 
from IKirsh et al.l (120091) is an estimate that is accurate only to within a factor 



of two (cf. IKirsh et al.l . 120091 Fig. 6), which is consistent with our results. 



We can also estimate the largest radius of mono-dispersed planetesimals that 
results in outward migration of the embryo, using our understanding of gas- free 
planetesimal - driven migration and aerodynamic drag. We recall that in the 
gas-free case, planetesimal scattering preferentially induces inward migration 
of the embryo. This arises from the slight imbalance of angular momentum 
transferred to the embryo by scattering planetesimals from either side of its or- 
bit, leaving the embry o with a dear t h of a ngular momentum and hence causing 
it to migrate inwards (IKirsh et al.l . 120091 ). With gas present, the aerodynamic 
gas drag will cause the orbits of planetesimals to circularize, thereby removing 
the planetesimals from the embryo's encounter zone. 



So we propose that the two relevant timescales are the timescale for the aero- 
dynamic gas drag to damp a planetesimal's eccentricity e by a factor of order 
itself (i.e. T eiaero = |e/e| aero ), and the migration timescale of the embryo due 
to planetesimal scattering (i.e. r a)Sca = \a/a\ sca ). Equating the two timescales 
determines a rough scaling relation for the transition between outward and 
inward embryo migration. 



If we substitute Eq. 10 for the gas volume density for a MMSN in the mid-plane 
of the disk and the Keplerian velocity, the aerodynamic gas drag timescale 
Eq. 4 can be written: 



where we have substituted in for the values of p p i sm i = 0.5 g cm -3 and Cd 
0.5. 

Using Eq. 5a, Eq. 18 and assuming / = e/2, we can write: 



2.3 x 10 2 year /l.0\ (R A ^\ ( a \ 13 / 4 



(19) 
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where f(a,e) = ^(a/AU) + 7.5 x 10- 5 (M em /M e ) 2 / 3 e| with r] = 1.95 x 
1CT 3 AlP 1 / 2 and = e/£/, is the Hill eccentricity. 

Using Eq. 2, we can define a migration timescale due to planetesimal scatter- 
ing: 



r a ,sca * 4.8 x 10 4 year ( M j ( JL) (20) 

Comparing the timescales given in Eq. 19 and Eq. 20, we define the maximum 
planetesimal radius (i.e. -Rdrag < -Rdra g ,trans) that results in outward migration 
of an embryo through the criterion r eiacro < r ajSca : 




-9/4 

iWrans " 2 -l X 10* km f(a, e) I ^ ] ( ) (21) 



In Fig. 7 we have plotted as a function of a the values of -Rdra g ,trans from 
Eq. 21, for four different e rms . Also plotted are the results from a series of 
simulations that were produced to determine the transition location. We see 
that this simple scaling relation is able to predict the transition to outward 
embryo migration. 

In the next section we create a more comprehensive toy model that addresses 
the different migration regimes discussed in this section: inward migration, 
outward migration and gap clearing. The model will make predictions for the 
transitions between these different migration regimes. 



5 Toy Model 



In this section we consider a toy model which displays the main features of 
the dynamical regimes described above. Indeed, with representative choices 
of the parameters it gives rough quantitative predictions for the planetesimal 
size for which outward migration is most likely to be triggered for a single 
embryo embedded in a dynamically cold planetesimal disk when aerodynamic 
drag effects are included. 



We begin by recalling some key results described in iKirsh et al.l (120091 ) for the 
gas-free case: 

(1) As noted in section §3.1, for sufficiently low mass embryos, planetesimal 
- driven migration is self-sustaining in either direction at a rate given by 
Eq. 2. 
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Fig. 7. Comparison of analytical estimate of the maximum iZdrag that results in 
outward migration from Eq. 21 in a fiducial MMSN disk, with results from simula- 
tions. Plotted are curves of -Rdrag,trans for different assumed RMS eccentricity e rms 
of the planetesimal population, in units of the Hill eccentricity e^. The points and 
error bars are estimates of -Rdrag,trans , extrapolated from five sets of 16 identical 
simulations at each location in the disk. The value of the drag coefficient is fixed at 
Cd = 0.5 in these simulations. 

(2) The steady-state migration develops from an instability caused by a pos- 
itive feedback loop in which the migra tion initially accelerates at a rate 



proportional to the migration rate (cf. iMasset and Papaloizoul . 120031 . for 
analogous result in Type - III migration in gas disks). 

(3) Unless there is a very strong positive disk density gradient in the outward 
direction, inward migration is the typical outcome. 

(4) The trigger for the inward migration arises from a slight asymmetry be- 
tween the rates that angular momentum is transferred to the embryo by 
outward versus inward scattering of planetesimals. 

The last point is demonstrated using results from the circular restricted three- 
body problem, in which the planetesimals have negligible total mass com- 
pared to th e embr yo. Starting with an initial dynamically cold population, 



Kirsh et al.l (120091 ) monitored the populations in the so-called "encounter" 
zones on either side of the embryo, which are the two regions of semi-major 
axes extending from roughly 1.0 - 3.5 Rh interior and exterior to the planet. 
Because of the existence of a constant of the motion known as the Jacobi in- 
tegral, only particles in the encounter zones can come within the Hill sphere 
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of the planet and have a strong scattering (or be accreted) . In iKirsh et al. 



( 120091 ) it is shown that there is an asymmetry in the transfer rates between 
the two reservoirs that leads to a build-up of particles exterior to the embryo. 
When the planetesimals have non-zero mass, the build-up of exterior particles 
corresponds to an outward flux of planetesimal angular momentum which, by 
conservation of angular momentum, causes the embryo to move inward. As 
fresh planetesimals interior to the embryo enter the feeding zone, the transfer 
asymmetry removes more angular momentum from the embryo, triggering a 
positive feedback i.e. an instability. Although the entire process is a compli- 
cated series of scatterings together with a slow diffusion in the planetesimal 
eccentricities, our simple model is parameterized by a single transfer timescale 
for each reservoir, together with a single estimate of the specific angular mo- 
mentum exchanged with the embryo as a particle is transferred from one side 
to the other. 

Thus our toy model has two reservoirs - one interior to the planet containing 
total time dependent mass M- m t(t) and one exterior with total mass M ext (t). 
We define AJ to be the gain of specific angular momentum of a particle 
scattered from the interior to the exterior reservoir. In the model we assume 
the planet of mass M em is fixed on a circular orbit of semi-major axis a em and 
infer the gain/loss in its angular momentum J em by equating it to the negative 
of the loss/gain of angular momentum of the scattering planetesimals. In the 
gas free case, then, the relevant equations for the toy model are: 
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(22a) 
(22b) 
(22c) 



where r in t and r cxt are the timescales for a planetesimal to be scattered out of 
the interior or exterior reservoir, respectively. 

The choice of the parameters in the m odel may be deter mined from restricted 
three-body scattering experiments (cf. IKirsh et all 120091 Fig. 1 and Fig. 2 and 



accompanying discussion). The initial asymmetry in the inner and outer en- 
counter zones means M ext (0) = M int (0)(l-|-5.5^/ l ). The timescale for scattering 
a particle from the inner zone is the product of the timescale for each encounter 
(i.e. the synodic period) time s the probab i lity p er encounter that such a scat- 
tering occurs. From Fig. 1 of IKirsh et al.l ( 120091 ) we find r in t = 6.7/^ orbital 
periods to be representative of a particle which starts near the middle of the 



encou nter zone. The value of r ext can be inferred from Fig. 2 of IKirsh et al. 



( 120091 ). together with the fact that the asymptotic solution of the equations 
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above is r cxt /r in t = M ex t/M int : we find r ext /r int = (1 + 12^). 

Finally, the value of A J can be estimated by noting that the conservation of 
the Jacobi constant during an outward scattering event typically transforms a 
particle with eccentricity e and apocentre near a em — Rh to one with pericentre 
near a em + Rh and comparable eccentricity. It is easy to show that for small 
eccentricity, this results in a change of specific angular momentum equal to 
(61 + e) J em , where J em = ^GM Q a em is the specific angular momentum of a 
circular orbit at the planet's location. 



The ecce ntricities of p articles in the encounter zones are typically a few £h- 
However, iKirshl ( 20071 ) found that particles within a given encounter zone dif- 
fusing to larger eccentricities (due to successive small scatterings) also trans- 
ferred angular momentum to the embryo. He found that the particles in the 
exterior reservoir evolved to higher eccentricities than those in the interior 
reservoir. Thus the net effect of diffusion in the two reservoirs results in a 
net gain in angular momentum of the planetesimals. Indeed he found the net 
angular momentum transferred by particles jumping from one reservoir to the 
other to be comparable in magnitude and sign to the net angular momen- 
tum transferred by diffusion within the reservoirs. This is incorporated in the 
model by using a somewhat larger value of the angular transfer rate than 
would be the case for particles that are only scattered between wings. We find 
a representative choice to be A J = 8£,hJem- 



Fig. 8 shows the behaviour of our toy model for the case of a 1.0M e planet 
on a circular orbit at 5.0 AU embedded in a sea of massless planetesimals. 
The upper solid curve shows the ratio predicted by Eq. 22 for the number 
of particles in the outer zone versus those in the inner encounter zone. The 
bottom solid curve shows the predicted average amount of angular momentum 
transferred per particle in units of the specific circular angular momentum. In 
both plots the results of direct iV-body simulations are shown for comparison. 
It can be seen that the agreement is very good, as it is for several other cases 
which we ran with different mass planets. 

Let us now consider the net angular momentum transferred to the embryo 
when aerodynamic drag acts on the planetesimals. Since the aerodynamic 
drag tends to circularize orbits, planetesimals scattered outward onto eccen- 
tric orbits tend to have their perihelion distance drawn from near the Hill 
sphere of the embryo to locations further away from the strong scattering re- 
gion. Indeed if the eccentricity damping timescale is shorter than the timescale 
for another strong scattering, the particle may evolve into a near-circular orbit 
several Hill radii from the embryo where it can be trapped in a mean-motion 
resonance exterior to the embryo and no longer have strong scatterings. On 
the other hand, planetesimals scattered inward will have both their aphelion 
distances and semi-major axes drawn inward. For rapid damping, planetes- 
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Fig. 8. Upper panel: the evolution with time of the ratio of masses in the two 
reservoirs on either side of a 1.0 M® planet at 5.0 AU. Lower panel: the average 
change in angular momentum of the planetesimals versus time. The solid curves 
are the prediction of the toy model discussed in the text. The dashed curves are 
the results of eight iV-body simulations of the same system using different random 
seeds for the initial planetesimal positions. 

imals scattered inward may then be dragged inward and not have another 
strong scattering with the embryo. The aerodynamic drag can be thought 
of as producing a "sink" for planetesimals from the two reservoirs. Moreover 
since the aerodynamic drag depends on the gas density, which itself has a 
strong radial dependence, there will be significant differences in the timescales 
for the draining of planetesimals from the interior versus exterior reservoir. 
In our simple toy model we thus add two terms associated with the removal 
of particles due to aerodynamic drag and drag timescales Td ra g,mt and Td ra g,cxt, 
for each reservoir respectively. The equations for our toy model now become: 



M int = M sca - 



M 
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M. 



ext 



7drag,int 

M ext 



7~drag,ext 



Jem M^sca^J 



(23a) 

(23b) 
(23c) 



where M sca = — M- int /T int + M cxt /r cxt is the net mass loss rate due to planetes- 
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imal scattering for the interior reservoir. 



Although analytic solutions to Eq. 23 exist, for the purposes of computing 
specific angular momentum transferred to the planet simple numerical inte- 
grations illustrate the qualitative nature of the solutions, and these can indeed 
be quantitatively compared to the results of the full iV-body simulations. For 
illustrative purposes, we adopt the gas disk density profile of LTD 10, in which 
the gas volume density is 3.4 x 10~ 9 g cm -3 at 1.0 AU and a = 9/4. For 
simplicity, in both the toy model and the iV-body simulations we fix the 
drag coefficient Co = 0.5. We insert a planet with mass M cm = 1.0 M e 
on a circular orbit at 5.0 AU in a planetesimal disk with surface density 
such that the initial mass in planetesimals in the inner encounter zone is 
-^int(O) = 2.0 M em . A typical particle in the scattered wings will have e ~ 3^ 
and using the definitions above we find for a planetesimal with drag radius 
i?dra g that T dragi i n t = 674 year (.R drag /km) and r dragiCxt = 1382 year (i? drag /km). 

Experience with the iV-body simulations shows that once the planet has mi- 
grated a distance ~ -R/i/3 in one direction the instability has been initiated 
and the planet continues to migrate in that direction. Thus we define: 



as the critical amount of angular momentum required to initiate self-sustaining 
migration. 

Plotted in Fig. 9 are the predictions of the toy model for the time evolution 
of the angular momentum transferred to the planet for the disk parameters 
given above and for several values of the planetesimal drag radius. For each 
value of -Rdrag the curve is terminated if the magnitude of the transfer equals 
A J crit since at that time self-sustaining migration (either inward or outward) is 
assumed to ensue. An alternative way to view the predictions of the toy model 
is shown in Fig. 10, in which the first extremum in the angular momentum 
transferred is shown as a function of -Rdrag- There it can be seen that for 
-Rdrag < 0.6 km a gap is predicted because insufficient angular momentum is 
transferred to trigger migration, and since in that case the encounter zones 
become evacuated, a gap in the planetesimal disk forms around the planet. 
For -Rdrag > 20 km, sufficiently negative angular momentum is transferred in 
the early stages to trigger inward migration. For values in between these drag 
radii outward migration is predicted. The discontinuity near -Rdrag = 20 km 
corresponds to a transition between solutions which "bottom-out" at A J em < 
— AJ crit and those which do not. 

In full iV-body simulations, individual scattering is of course stochastic and 
the fairly subtle asymmetries involved therefore lead to stochastic outcomes. 




cm 




(24) 
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Fig. 9. The time evolution of the angular momentum transferred to the planet for 
the toy model parameters discussed in the text. Thicker curves correspond to larger 
planetesimals, as shown in the legend. 

To compare with the toy model for the same disk parameters we have run 
numerous iV-body experiments and binned the outcomes into 3 possible sce- 
narios: outward migration, inward migration and planetesimal gap formation. 
The simulations involve a 1.0 M e embryo at 5.0 AU, with a population of 
3.0 x 10 4 dynamically cold (i.e. a e = 2a i = 10~ 3 ) planetesimals distributed 
from 4.0 AU to 6.0 AU, each with a different Rd rag value. Approximately 16 
A^-body simulations were performed, in each of 10 logarithmically spaced bins 
in .Rdrag- Each simulation involved a 1.0 M e embryo at 5.0 AU, embedded in a 
sea of massless planetesimals. We note that Fig. 10 and Fig. 11 agree with each 
other, at least on a qualitative level. While the toy model is able to capture 
the physics underlying embryo migration, it fails to account for the minutia 
of the planetesimals' dynamical behaviour. So while our toy model can rea- 
sonably predict the planetesimal radius which most often leads to outward 
embryo migration, it also predicts values for i?drag,tra P and i?dra g ,trans that are 
almost an order of magnitude too small and too large, respectively. 

We can extrapolate the results in Fig. 11 as a function of M em and a em via 
the transitional radii denned by Eq. 16 and Eq. 21. If we assume Cd = 0.5 
for simplicity then for fixed M em , i?dra g ,tra P and -Rdra g ,trans will both decrease 
with increasing a em . However, since i?dra g ,tra P decreases faster than -Rdrag.trans, 
the width of the regime for outward embryo migration will increase with a em . 
In the case for fixed a em , i?dra g ,tra P will decrease, while i?dra g ,trans will increase 
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with M em . So the width of the regime for outward embryo migration will 
expand with M em , albeit differentially. While an extrapolation of the results 
in Fig. 11 could also be performed with aC^ computed based on the dynamical 
and physical properties of a typical planetesimal, the nonlinear nature of Co 
does not permit for simple scaling relations. In addition, extensive numerical 
simulations to confirm these predictions would be computationally expensive. 

To demonstrate the effect of changing the slope of the gas density profile, 
Fig. 12 shows the results of iV-body simulations of a 1.0 M e mass planet at 
5.0 AU in which the gas volume density is 9.0 x 10~ n g cm -3 at 5.0 AU, 
Cd = 0.5 and r\ is determined as in §3.2.1. By varying a and keeping the gas 
to solid ratio fixed {e.g. (5 = a — 5/4) we can see the important effects on 
outward migration that a large gradient in gas density creates. Recall that the 
typical semi-major axes of particles in the exterior reservoir, are somewhat 
larger than in the interior reservoir. Since the gas drag timescale is inversely 
proportional to the local gas density, a strong gradient in gas density thus leads 
to a large disparity in the gas drag timescale between the two reservoirs. Thus 
particles in the inner reservoir are differentially depleted more rapidly when 
the gas density gradient is large. The results of the toy model (not shown) are 
in good qualitative agreement with the iV-body simulations. 
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Fig. 11. Summary of results from multiple iV-body simulations of the evolution 
of a 1.0 M® embryo starting at 5.0 AU for the disk parameters described in the 
text. Three curves are shown: for each value of planetesimal drag radius -Rdrag the 
solid curve shows the fraction migrating out, the dot-dash curve shows the fraction 
migrating in and the dashed curve displays the fraction migrating in. 

6 Size Distribution 



Up to this point we have considered only mono-dispersed planetesimals in our 
simulations, but as was noted earlier this is not very realistic. While incor- 
porating a self-consistent time dependent planetesimal size distribution is not 
feasible at present, we did want to address concerns about the combined ef- 
fect of a size distribution on embryo migration. To that end, we implement a 
static planetesimal size distribution: dN/ dR^ ra _ g oc -R^i-ag for planetesimal sizes 
0.0625 km < -Rdrag < 32.0 km, and we investigate the behaviour of the embryo 
migration for a range of power-law exponents: 2 < 7 < 5. Of special note are 
two values of the power-law exponent, namely 7 = 7/2 and 7 = 4. In the case 
of 7 = 7/2, this value correspond s to the size distribution resulting from a col- 



lisional cascade (IDohnanyil . Il969[ ) . while the value of 7 = 4 corresponds to the 



case where the mass in the size distribution is logarithmically equally-spaced. 

It stands to reason that if a significant fraction of the mass in the planetesi- 
mal size distribution preferentially causes the embryo to migrate in a certain 
direction, then that population of planetesimals will dictate the outcome of 
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Fig. 12. Similar to Figure 11 but for iV-body simulations of a 1.0 M© mass planet at 
a = 5.0 AU in which the gas density is 9.0 x 10~ n g cm~ 3 (a/5.0 AU)~ a , Co = 0.5, 
and 77(a) = 6.0 x 10~ 4 (a + ^(a/AU) 1 / 2 . From top to bottom, the results shown are 
for a = {0.0, 1.0, 2.25, 3.25}, respectively. 

embryo migration. We can estimate the range of 7 that will give rise to the 
different migration outcomes (cf. §4.1) by computing the fraction of the total 
mass those -Rdrag will contribute that give rise to the migration outcomes (cf. 
Fig. 13). For an embryo of mass 0.25 M e at 10.0 AU in a fiducial MMSN 
disk, the gap clearing regime occurs for -Rdrag < 0.125 km. For the case of 
outward migration, that occurs for 0.125 km < -Rdrag ^ 2.0 km, and finally for 
inward migration occurs for -Rdrag > 2.0 km. The fraction of the total mass 
that these different regimes contribute are plotted in Fig. 13. What this plot 
indicates is that there will be a smooth transition between the different migra- 
tion regimes, but since the migration outcome is bimodal, stochastic effects 
will be important. 

To determine what occurs when we include a size distribution, we perform 16 
sets of simulations for a 0.25 M e embryo at 10.0 AU in a fiducial MMSN disk, 
each for seven values of 7. The average migration rate for the embryo from 
these simulations are plotted in Fig. 14 as a function of 7. We see that for 
4.0 < 7 < 4.5, the average embryo migration rate is outwards while for all 
other values of 7 the embryo migrated inwards. In Fig. 15 we plot the fraction 
of the 16 simulations that result in outward migration, and we note that for 
7 < 3.5 there are still ~ 20-40% of the simulations that result in outward 
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Fig. 13. The fraction of the total mass contained in planetesimals that lead to one 
of the migration outcomes: inward migration, outward embryo migration or gap 
clearing, as a function of the exponent 7 for the planetesimal size distribution. For 
a 0.25 M(q embryo at 10.0 AU in a fiducial MMSN disk, the transitions between 
these three different outcomes occurs at: ~ 0.125 km and ~ 2.0 km. This assumes a 
range 0.0625 km < i?drag < 32.0 km for the planetesimal size distribution. 



migration. This is despite the fact that these planetesimal size distributions 
heavily favour large Rdr&g, and consequently inward migration. If we now define 
the 50% fraction as the threshold defining outward embryo migration, then 
the range of 7 that results in outward migration is: 3.7 < 7 ou t ^ 4.7. For this 
range of 7, the majority of the mass in size spectrum are in those planetesimals 
that cause the embryo to migrate outward. So for disk models with a different 
range of planetesimal radii causing outward migration of the embryo, then 
we would expect the range of 7 to shift to smaller or larger values depending 
if the planetesimal radii "sweet spot" makes the size spectrum less or more 
top-heavy. 

So we can see that for a reasonable range of 7 the planetesimal size distribution 
results in outward embryo migration, but even outside this range there is 
still a small chance that the embryo will migrate outwards. Next we look at 
extending our analysis to include Type - I migration, and see how it competes 
with planetesimal - driven migration and aerodynamic gas drag. 
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Fig. 14. Migration rate for an 0.25 M® embryo at 10.0 AU in a fiducial MMSN disk, 
as function of the exponent 7 of the planetesimal size distribution. The average of 16 
identical simulations were used to determine the migration rate and its uncertainty, 
where only the initial planetesimal distribution were randomized. The value of Cd 
is computed based on the dynamical and physical properties of each planetesimal in 
these simulations, and a range of planetesimal radii: 0.0625 km < -Rdrag < 32.0 km 
were used. 

7 Effects of Type - I 



We now turn our attention to Type - I migration, and address its ability to 
deposit embryos and giant planetary cores on the central star. We could ar- 
gue that this tendency implies that the efficiency parameter c a as it appears in 
Eq. 8, must be much less than unity for large embryo masses. While decreasing 
c a would solve the issue with embryo deposition, it cannot be justified in gen- 
eral except where the physical conditions permit it (e.g. depleted gaseous disks, 
turbulence, etc.). As we have shown in §4.1, planetesimal - driven migration 
can lead to rapid inward and outward embryo migration. Furthermore, simu- 
lations with outward migration also result in the rapid growth of the embryo, 
particularly in dynamically cold disks. An illustration is provided in Fig. 16, 
where a 0.25 M® embryo is placed at 5.0 AU in disk similar to the example 
used in §5 with f g — f s — \/5 and a = 9/4 (LTD10) with a planetesimal size 
-Rdrag — 1.0 km. At the end of the simulation, the embryo migrates outwards 
to ~ 14.5 AU and grows to a mass of ~ 8.5 M®. This example suggests that 
including a population of cold, collision fragments, especially if planetary en- 
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Fig. 15. Fraction of simulations that resulted in outward migration as a function 
of the exponent 7 of the planetesimal size distribution, for a 0.25 M ffi embryo at 
10.0 AU in a fiducial MMSN disk. Sets of 16 identical simulations were used to 
determine the outward fraction, where only the initial planetesimal distribution 
were randomized. In these simulations. The value of Co was computed based on 
the dynamical and physical properties of each planetesimal in these simulations, 
and a range of planetesimal radii: 0.0625 km < -Rdrag < 32.0 km were used. The 
/out = 0.5 threshold (dashed line) was used to determine range of 7 that results in 
outward migration. 



velopes are included (jlnaba and Ikomal . 120031 ) could aid in the rapid formation 
of a core. 



Under conditions for outward planetesimal - driven embryo migration and 
growth, Type - I migration will act to counteract this outward migration. It 
stands to reason with the embryo accreting planetesimals, there will be a point 
where Type - I migration will dominate planetesimal - driven migration. To 
assess the relative importance of these two competing migration mechanisms, 
we start by asking when would the timescale for Type - I migration T a)typc i 
be shorter than the timescale for planetesimal - driven migration r a sca . The 
formula for these two timescales are given below, and in the case of T a t y pcI we 
assume e < 2^ h . 



7~a,typel 



Povh ( %s\ 



2-KC a \ a 



S gas (a)7ra 2 

M 
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Fig. 16. Simulation of a 0.25 M® embryo at 5.0 AU in a disk with f g = f s = \/5, 
a = 9/4 and 1.0 km size planetesimals. After 10 5 yrs, the embryo migrated outwards 
to ~ 14.5 AU and attained a mass of ~ 8.5M ffi . 
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We want to know at which embryo mass M em will Type - I dominate the 
planetesimal - driven migration. Equating Eq. 25 and Eq. 26, and solving 
the resulting cubic equation for M em , we find that this occurs when M em > 

-^em,crit • 



G{a) 

W) 

where the functions F(a) and G(a) are given below: 



M em>crit = AF(a) - B 



(27) 



F(a) 
G(a) 
H(a) 



dH(a) + C 2 y/C 3 H 2 (a) + C A G 3 (a) 



1/3 



1.0 / 
1.0 



5.0 AuJ 



c a 
1.0 



60.0 
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\ 0.05/ 
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where f gs is the gas to solid ratio, ho is the aspect ratio of the disk at 1.0 AU, 
and the coefficients are listed in Table 1. 



A (M e ) 


B (M e ) 


Ci 


c 2 


c 3 


c 4 


0.3816 


0.2522 


26.50 


1.732 


2.341 x 10 2 


9.623 x 10"2 



Table 1 

Coefficients for Eq. 27 and Eq. 28. 



To test the predicted values of M em ci . it , we performed 8 sets of simulations of 
an embryo migrating outwards for several different f s and c a values, accreting 
planetesimals as it migrates. We assess the critical mass to be when the embryo 
ceases to migrate outward, and reverses direction. The simulations start with 
a 0.25 M e embryo at 5.0 AU, with a population of 6.4 x 10 4 dynamically cold 
(i.e. a e = 2a i = 10 -3 ) planetesimals distributed from 5.0 AU to 15.0 AU. All 
these simulations are integrated for 10 5 year with a 0.5 year time step, and we 
assume the Type - I efficiency parameter for eccentricity is set to c e = 1.0. 
Plotted in Fig. 17 is the maximum embryo mass attained for disk models with 
f t = {0.5, 1.0, 2.0}^ and a = 9/4, each for c a = {0.1, 0.5}. 

First we note that the typical maximum mass attained is below the desired 
10.0 Me, however the range of maximum mass does approach this limit for 
larger f s values, and with marginal dependence on c a . We also note that the 
maximum mass attained in these simulations do not exceed the predicted 
values from Eq. 27, even within the error bars. This may be a consequence of 
the relatively low efficiency of planetesimal accretion, or that the criterion does 
not capture all the dynamics. However in LTD10 they demonstrated that the 
addition of an atmosphere can dramatically increases the accretion efficiency, 
which is a promising prospect. Though this atmosphere code can be made 
available, we shall leave the inclusion of atmospheres for future work. 



We have omitted the c a = 1.0 case in this investigation since it can be demon- 
strated that for a viscously evolving disk, one would need to consider an em- 
bryo mass an order of magnitude smaller than inves tigated in this stu dy. This 



can be demonstrated by using the results found in iLyra et al.l (120101 ). where 



simulations of viscously and radiatively evolving protoplanetary disks are pre- 
sented. From these simulations (See their Fig. 3) we assess the viscous accre- 
tion timescale of their disks to be on the order of ~ 5 Myr, and if this timescale 
is compared to the Type - I timescale in Eq. 25, we can place a constraint 
on M em or c a . Assuming the fiducial values for the disk model we consider in 
our simulations, this gives M em ~ 0.02 M @ /c a . This relation implies that we 
would need to consider a 0.02 M e embryo if c a = 1.0, which is an embryo mass 
an order of magnitude smaller than we consider in this work. 



As a final note, while the eventual dominance of Type - I migration does not 
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Fig. 17. The critical embryo mass as a function of the solid surface density scaling 
factor f s , above which the Type - I migration timescale is shorter than the plan- 
etesimal - driven migration timescale (cf. Eq. 27). This assumes f g = f s for a disk 
model with a = 9/4 at 5.0 AU. The curves correspond to M emiCr i t for three different 
values of of the Type - I efficiency factor c a , while the points are the average maxi- 
mum mass attained from simulations. The simulations start with a 0.25 embryo 
at 5.0 AU, and where the circles and squares are for simulations with c a = 0.1 and 
c a = 0.5, respectively. 



solve the problem of cores being deposited onto the central star, it is important 
to note the location of the giant planetary core when Type - I migration 
begins to dominate. Since planetesimal - driven migration can migrate an 
embryo several AU further away from the central star, this can significantly 
increase the timescale for the giant planetary core to reach the central star. 
For the example illustrated in Fig. 16, tripling the embryo's initial distance 
to ~ 14.5 AU increases the Type - I migration timescale in Eq. 25 by more 
than an order of magnitude (e.g. T aitypc i ~ 1.4 Myr). This may be sufficient 
to permit the giant planetary core to survive long enough for it to accrete a 
massive gaseous envelope, and with the dispersal of the gaseous disk decrease 
the effectiveness of Type - I migration to deposit planets onto the central star. 
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8 Conclusions 



To explain the observed distribution of giant extrasolar planets, including the 
giant planets in the Solar System, a comprehensive model of their formation 
is needed. However, the physical processes and conditions involved give rise to 
very complex dynamical behaviour. In view of the results of LTD10, a detailed 
study of planetesimal scattering and gas drag fills a much needed niche. The 
resulting dynamics for a single embryo embedded in a sea of planetesimals 
subject to aerodynamic drag is quite fascinating and provides clues to future 
avenues of research. 

The findings of this research are: 

(1) The presence of a gas disk modifies planetesimal - driven migration of 
a single embryo, resulting in four different migration regimes: streaming, 
gap clearing, outward and inward embryo migration. 

(2) The range of planetesimal radii that result in outward embryo migration 
are plausible, and represents a possible dynamical avenue for embryo 
growth. 

(3) There exists simple scaling relations for the transition between the stream- 
ing and gap clearing regimes, along with the outward and inward embryo 
migration regimes. 

(4) Our toy model provides a good qualitative understanding of the observed 
phenomena in the A-body simulations. 

(5) Simulations with a static planetesimal size distribution also demonstrate 
the different migration outcomes, and it was shown that outward migra- 
tion is possible for a finite range of the power-law exponent of the size 
distribution. 

(6) Inclusion of Type - I migration with planetesimal - driven migration and 
aerodynamic gas drag were explored. An analytical estimate was shown of 
the maximum achievable embryo mass before Type - I migration begins 
to dominate over outward planetesimal - driven migration, which appears 
to be borne out in simulations. 

As described in a companion paper (LTD 10), future work will incorporate 
other effects: multiple embryos, evolving planetesimal size distribution and 
embryo atmospheres, and the effects of fragmentation. 



Acknowledgements 

The ongoing financial support of NSERC is gratefully acknowledged. HFL is 
grateful for funding from NASA's Origins and OPR programs, as well as a 



34 



grant from the National Science Foundation (Award ID 0708775). 



References 

I. Adachi, C. Hayashi, and K. Nakazawa. The gas drag effect on the elliptical 

motion of a solid body in the primordial solar nebula. Progress of Theoretical 

Physics, 56:1756-1771, 1976. 
R. Brasser, M. J. Duncan, and H. F. Levison. Embedded star clusters and the 

formation of the Oort cloud. II. The effect of the primordial solar nebula. 

Icarus, 191:413-433, 2007. 
E. I. Chiang and P. Goldreich. Spectral Energy Distributions of T Tauri Stars 

with Passive Circumstellar Disks. ApJ, 490:368-+, 1997. 
J. N. Cuzzi and K. J. Zahnle. Material Enhancement in Protoplanetary Neb- 
ulae by Particle Drift through Evaporation Fronts. ApJ, 614:490-496, 2004. 
J. S. Dohnanyi. Collisional Model of Asteroids and Their Debris. Journal of 

Geophysical Research, 74:2531 — h, 1969. 
M. J. Duncan, H. F. Levison, and M. H. Lee. A Multiple Time Step Symplectic 

Algorithm for Integrating Close Encounters. A J, 116:2067-2077, 1998. 
J. A. Fernandez and W.-H. Ip. Some dynamical aspects of the accretion of 

Uranus and Neptune - The exchange of orbital angular momentum with 

planetesimals. Icarus, 58:109-120, 1984. 
P. Goldreich and S. Tremaine. Disk-satellite interactions. ApJ, 241:425-441, 

1980. 

R. S. Gomes, A. Morbidelli, and H. F. Levison. Planetary migration in a 
planetesimal disk: why did Neptune stop at 30 AU? Icarus, 170:492-507, 
2004. 

K. E. Haisch, Jr., E. A. Lada, and C. J. Lada. Disk Frequencies and Lifetimes 
in Young Clusters. ApJL, 553T153-L156, 2001. 

C. Hayashi. Structure of the Solar Nebula, Growth and Decay of Magnetic 
Fields and Effects of Magnetic and Turbulent Viscosities on the Nebula. 
Progress of Theoretical Physics Supplement, 70:35-53, 1981. 

O. Hubickyj, P. Bodenheimer, and J. J. Lissauer. Accretion of the gaseous 
envelope of Jupiter around a 5 10 Earth-mass core. Icarus, 179:415-431, 
2005. 

S. Ida and J. Makino. N-body simulation of gravitational interaction between 
planetesimals and a protoplanet. I - Velocity distribution of planetesimals. 
Icarus, 96:107-120, 1992. 

S. Inaba and M. Ikoma. Enhanced collisional growth of a protoplanet that 
has an atmosphere. A&A, 410:711-723, 2003. 

D. M. Kary, J. J. Lissauer, and Y. Greenzweig. Nebular gas drag and planetary 
accretion. Icarus, 106:288 — h, 1993. 

D. R. Kirsh. Simulations of planet migration driven by the scattering of smaller 
bodies. Master's thesis, Queen's University (Canada), 2007. 



35 



D. R. Kirsh, M. Duncan, R. Brasser, and H. F. Levison. Simulations of planet 
migration driven by planetesimal scattering. Icarus, 199:197-209, 2009. 

E. Kokubo and S. Ida. Oligarchic Growth of Protoplanets. Icarus, 131:171- 
178, 1998. 

D. G. Korycansky and J. B. Pollack. Numerical calculations of the linear 
response of a gaseous disk to a protoplanet. Icarus, 102:150-165, 1993. 

H. F. Levison, E. Thommes, and M. J. Duncan. Modeling the Formation of 
Giant Planet Cores. I. Evaluating Key Processes. A J, 139:1297-1314, 2010. 

W. Lyra, S.-J. Paardekooper, and M.-M. Mac Low. Orbital Migration of 
Low-mass Planets in Evolutionary Radiative Models: Avoiding Catastrophic 
Infall. ApJL, 715:L68-L73, 2010. 

R. Malhotra. The origin of Pluto's peculiar orbit. Nature, 365:819-821, 1993. 

F. S. Masset and J. C. B. Papaloizou. Runaway Migration and the Formation 
of Hot Jupiters. ApJ, 588:494-508, 2003. 

H. Mizuno. Formation of the Giant Planets. Progress of Theoretical Physics, 
64:544-557, 1980. 

S.-J. Paardekooper and J. C. B. Papaloizou. On disc protoplanet interactions 
in a non-barotropic disc with thermal diffusion. A&A, 485:877-895, 2008. 

S.-J. Paardekooper and J. C. B. Papaloizou. On corotation torques, horse- 
shoe drag and the possibility of sustained stalled or outward protoplanetary 
migration. MNRAS, 394:2283-2296, 2009. 

S.-J. Paardekooper, C. Baruteau, A. Crida, and W. Kley. A torque formula 
for non-isothermal type I planetary migration - I. Unsaturated horseshoe 
drag. MNRAS, 401:1950-1964, 2010. 

J. C. B. Papaloizou and J. D. Larwood. On the orbital evolution and growth 
of protoplanets embedded in a gaseous disc. MNRAS, 315:823-833, 2000. 

J. B. Pollack, O. Hubickyj, P. Bodenheimer, J. J. Lissauer, M. Podolak, and 
Y. Greenzweig. Formation of the Giant Planets by Concurrent Accretion of 
Solids and Gas. Icarus, 124:62-85, 1996. 

R. R. Rafikov. Fast Accretion of Small Planetesimals by Protoplanetary Cores. 
A J, 128:1348-1363, 2004. 

S. G. Sousa, N. C. Santos, M. Mayor, S. Udry, L. Casagrande, G. Israelian, 
F. Pepe, D. Queloz, and M. J. P. F. G. Monteiro. Spectroscopic parameters 
for 451 stars in the HARPS GTO planet search program. Stellar [Fe/H] and 
the frequency of exo-Neptunes. A&A, 487:373-381, 2008. 

E. W. Thommes, M. J. Duncan, and H. F. Levison. Oligarchic growth of giant 
planets. Icarus, 161:431-455, 2003. 

S. Udry, D. Fischer, and D. Queloz. A Decade of Radial- Velocity Discoveries 
in the Exoplanet Domain. Protostars and Planets V, pages 685-699, 2007. 

W. R. Ward. Density waves in the solar nebula - Differential Lindblad torque. 
Icarus, 67:164-180, 1986. 

W. R. Ward. Protoplanet Migration by Nebula Tides. Icarus, 126:261-281, 
1997. 

G. W. Wetherill and G. R. Stewart. Accumulation of a swarm of small plan- 
etesimals. Icarus, 77:330-357, 1989. 



36 



J. Wisdom and M. Holman. Symplectic maps for the n-body problem. AJ, 
102:1528-1538, 1991. 



37 



